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The dependence of protoplanet migration rates on coorbital torques 
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ABSTRACT 

We investigate the migration rates of high-mass protoplanets embedded in accretion discs via two and 
three-dimensional hydrodynamical simulations. The simulations follow the planet's radial motion and 
employ a nested-grid code that allows for high resolution close to the planet. We concentrate on the 
possible role of the coorbital torques in affecting migration rates. We analyse two cases: (a) a Jupiter-mass 
planet in a low-mass disc and (b) a Saturn-mass planet in a high-mass disc. The gap in case (a) is much 
cleaner than in case (6). Planet migration in case (&) is much more susceptible to coorbital torques than in 
case (a). We find that the coorbital torques in both cases do not depend sensitively on whether the planet 
is allowed to migrate through the disc or is held on a fixed orbit. We also examine the dependence of the 
planet's migration rate on the numerical resolution near the planet. For case (a), numerical convergence 
is relatively easy to obtain, even when including torques arising from deep within the planet's Hill sphere, 
since the gas mass contained within the Hill sphere is much less than the planet's mass. The migration rate 
in this case is numerically on order of the Type II migration rate and much smaller than the Type I rate, 
if the disc has 0.01 solar-masses inside 26 AU. Torques from within the Hill sphere provide a substantial 
opposing contribution to the migration rate. In case (&), the gas mass within the Hill sphere is larger 
than the planet's mass and convergence is more difficult to obtain. Torques arising from within the Hill 
sphere are strong, but nearly cancel. Any inaccuracies in the calculation of the torques introduced by grid 
discretization can introduce spurious torques. If the torques within the Hill sphere arc ignored, convergence 
is more easily achieved but the migration rate is artificially large. At our highest resolution, the migration 
rate for case (6) is much less than the Type I rate, but somewhat larger than the Type II rate. 

Subject headings: accretion, accretion discs — hydrodynamics — planetary systems: formation, protoplanetary 
discs 



1. Introduction 

When the first planetary systems were discovered, migra- 
tion provided the natural explanation for the existence of the 
so-called "Hot Jupiters" (Lin, Bodenheimer & Richardson 
1996). For this explanation to hold, migration time-scales 
should be no longer than disc life-times of several million 
years (e.g., Haisch, Lada & Lada 2001). However, the migra- 
tion and planet formation processes are inter-related. Clearly, 
there would be complications and possibly difficulties in un- 
derstanding planet formation by a process whose time-scale 
is long compared to the migration time-scale. 



1 To appear in the Monthly Notices of the Royal 
Astronomical Society". A version with high resolution 
figures is available at 

http : //www . astro . ex . ac . uk/people/gennaro/publicat ions/ . 



In the case of giant planet formation by the core accre- 
tion process (e.g., Bodenheimer & Pollack 1986; Wuchterl 
1991a), the formation time-scale of about 10 7 years (Pollack 
et al. 1996; Tajima & Nakagawa 1997) is rather long com- 
pared to migration time-scales of about 10 5 years (e.g., Lin & 
Papaloizou 1986; Ward 1997) for a Jupiter-mass planet and 
about 10 6 years for an Earth-mass core (Tanaka, Takeuchi 
& Ward 2002; D'Angelo, Kley & Henning 2003; Bate et al. 
2003). 

However, it has been found recently (Rice & Armitage 
2003; Alibert, Mordasini & Benz 2004) that the effects of ac- 
cretion and migration of a planetary core can significantly 
reduce the time needed by the core to reach the mass neces- 
sary for the nucleated instability to occur (Wuchterl 1991b; 
Magni & Coradini 2004). Furthermore, several recent studies 
have suggested that additional effects may be of importance 
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to migration. These include thermal effects of the disc mate- 
rial near a planet (Morohoshi & Tanaka 2003; Jang-Condell 
& Sasselov 2004), effects of radial opacity jumps in the disc 
(Menou & Goodman 2004), effects of vortices induced by a 
planet (Roller, Li & Lin 2003), effects of turbulent fluctu- 
ations (Nelson & Papaloizou 2004), and effects of coorbital 
material (Masset & Papaloizou 2003, hereafter MP03). In 
the current study, we consider effects of coorbital material, 
along the lines of MP03. 

Corotation torques arise in the coorbital region. In the 
absence of dissipation or other time-dependent effects, the 
corotation torque is zero in a smooth disc. The reason is that 
in a steady-state, fluid elements circulate in closed orbits. 
Over a libration time-scale, a fluid element will gain and lose 
torque, but the result is zero average torque. Formally, the 
corotation region in linear theory gives rise to a torque that 
depends on the gradient of the disc vortensity (e.g., Goldreich 
& Tremaine 1979). This torque is properly interpreted as 
an "unsaturated" or maximal torque that arises over time- 
scales less than a libration time-scale or when the effects of 
viscosity are sufficiently large in steady-state. A derivation 
that includes nonlinear feedback shows that the steady-state 
corotation torque is indeed zero for a fluid in a smooth inviscid 
disc (Balmforth & Korycansky 2001; Ogilvie & Lubow 2003, 
see also Masset 2001). But, even the unsaturated corotation 
torque for typical planet-disc systems is somewhat smaller in 
magnitude than the other (Lindblad) torques present (Tanaka 
et al. 2002). Furthermore, for typical disc parameters, this 
torque is likely saturated (reduced to a smaller value), since 
the effects of turbulent viscosity are not sufficiently strong, at 
least in an alpha model description. 

The above-described analyses did not take into account 
the effects of the radial migration of the planet. This motion 
may cause enough asymmetry in the corotational flow that a 
net torque occurs, which may lead to a "runaway" situation 
(MP03). That is, the migration of the planet might cause a 
corotational torque that enhances the migration rate, which 
in turn further promotes asymmetry and leads to a stronger 
torque, etc. Examples of such a runaway phenomenon were 
reported in simulations by MP03. The most favourable cir- 
cumstances for such a process are expected when a planet 
interacts with a massive disc in which there is not a clean 
gap. 

In addition to the classical corotational torques that arise 
from nearly librating orbits, coorbital torques can also arise 
within the Hill sphere of the planet. Material flows into 
this region and forms a circumplanetary disc with shocks 
(Lubow, Seibert & Artymowicz 1999; D'Angelo, Henning & 
Kley 2002). 

Our previous studies did not allow the planet to migrate 
during the course of the simulation and therefore could not 
have found such a runaway migration. Numerical resolution is 
a key issue because densities near a planet are relatively high 
and fractionally small density errors there can give rise to 
large spurious torques. Bate et al. (2003) found that torques 
near the planet may contribute somewhat (~ 20 per cent) to 
the migration rate. However, that study lacked the resolution 
to reliably determine such torques. 

In this paper, we investigate if the torques exerted on a 



high-mass planet by a disc depend significantly on whether 
the planet is kept on a fixed orbit or allowed to migrate. 
We also investigate the possible role of torques due to ma- 
terial within the Hill sphere. We do this by means of two- 
dimensional (2D) and three-dimensional (3D) high resolution 
hydrodynamical simulations. A key feature of the code is 
that it allows high resolution to be achieved by means of 
nested grids that encompass a region around the planet as 
it migrates. With this code, we are able to examine the con- 
tribution of the material inside the planet's Hill sphere to the 
total torque on the planet. 

The outline of the paper is as follows. In Section 2 the 
physical model is described. In Section 3 we present an 
overview of the numerical procedures employed in these com- 
putations. The results of the calculations are provided in 
Sections 4 and 5. In Section 6 we present a discussion of 
these results and our conclusions. 

2. Description of the physical model 

It is generally believed that the interaction between a cir- 
cumstellar disc and a Jupiter-sized object can be studied by 
means of a two-dimensional approximation (Kley, D'Angelo 
& Henning 2001; D'Angelo et al. 2003). However, while this is 
possibly true when considering interactions occurring at Lind- 
blad resonance locations (i.e., at distances from the planet 
larger than a disc scale-height, H), it is not yet clear whether 
or to what extent this remains a valid assumption when deal- 
ing with other interactions occurring at coorbital locations 
(Masset 2002). Therefore, in this investigation we considered 
both 2D and 3D disc models. 

In the 2D geometry we employed a cylindrical coordinate 
frame {O; r, <f), z}, with the disc confined in the plane 2 = 0, 
whereas in the 3D geometry we used a spherical polar coor- 
dinate frame {O; R, 9, </>}. The rotational axis of the disc is 
either parallel to the z-axis or to the polar direction, 8 = 0. 
Both reference frames have their origin, O, on the star and 
rotate about the disc axis with an angular velocity ft and 
an angular acceleration fl, being this last vector also parallel 
to the disc axis. The magnitudes of fl and ft are specified 
later in this section. For the sake of clarity we point out that, 
whenever the variable r is used in the context of spherical po- 
lar coordinates, it indicates the distance from the rotational 
axis r = R sin 6. 

2.1. Equations of motion for the disc 

The hydrodynamical equations describing the disc evolu- 
tion are usually written in the conservative form for the radial 
and angular momenta. This can be derived from the Navier- 
Stokes equations for the velocities (see, e.g., Mihalas & Weibel 
Mihalas 1999, Chapter 3) and the continuity equation. Since 
the 2D equations in cylindrical coordinates can be formally 
derived from the 3D equations in spherical polar coordinates, 
we explicitly write them only for the latter reference frame. 
Indicating with p the mass density, with u = (un,ug,u$) 
the fluid velocity, and with u>a — + Q the absolute angu- 
lar velocity of the fluid around the disc axis (ur = u$), the 
equations of motion for the disc in conservative form can be 
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written as 
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^+V.(^ti) = -g-p^+fl B infl^ I (4) 

where 

(Cfli^e.^) = p{u R ,u g R,ll> a R 2 sin 2 6») (5) 
are the radial and angular momentum densities. Equations 
in 2D cylindrical coordinates can be obtained from equa- 
tions (1), (2), (4), and (5) by replacing p with the surface 
density E, using the appropriate expression for the divergence 
operator, dropping all terms that contain the velocity ug, and 
setting 6 = tt/2. 

Note that ^ is the absolute azimuthal angular momentum 
(density) of the fluid rather than that relative to the rotating 
reference frame. This basically means that the non-inertial 
terms arising from the rotation of the reference frame (i.e., 
Coriolis and angular velocity accelerations) are incorporated 
in the left-hand side of equation (4). This choice assures a 
better numerical treatment of the associated conservation law 
(Kley 1998). 

We adopted a locally isothermal equation of state by set- 
ting p = c 2 p (or p = c 2 E in 2D) . The sound speed, c s , is 
equal to the disc aspect ratio, H/r, times the Keplerian ve- 
locity, vk- We used a constant disc aspect ratio throughout 
the disc, implying that the temperature distribution scales as 
the inverse of the distance from the disc axis. 

Since self-gravity is ignored, the gravitational potential, $, 
only includes contributions from the star, the planet, and the 
non-inertial forces arising from the motion of the frame origin, 
O. Indicating the position vector of a fluid element as x and 
that of the planet as x p , the disc gravitational potential reads 



$ = 



GM, 



G M v 



- X p \ 2 + £' 



+ 



GM„ 



(6) 



where M, is the stellar mass, M p is the planet mass, and 
e is a smoothing length (see the discussion in section 2.4). 
The third term on the right-hand side of equation (6) origi- 
nates from the fact that the origin of the coordinate frame is 
accelerated by the planet 2 . 

The viscosity force density, / = (fn,fe,ftj,) (or / = 
(f r ,f<p) in 2D), is written as / = V • S. It assumes a stan- 
dard viscous stress tensor, S, for a Newtonian fluid with a 
constant kinematic viscosity, u, and a zero bulk viscosity. Ex- 
plicit forms for the components of / can be found in Mihalas 
& Weibel Mihalas (1999, Chapter 3), for the 3D spherical 
polar coordinates case and in D'Angelo et al. (2002), for the 
2D cylindrical coordinates case. 



2 To be strict, an additional term should appear in equation (6) 
due to the force exerted by the disc material on to the star, as 
measured from the centre-of-mass reference frame. We neglected 
this contribution, as is done when assuming that the centre-of-mass 
of the whole system coincides with that of the star— planet system. 



2.2. Equation of motion for the planet 

In the present study the planet's orbit evolves under the 
gravitational action of the central star and of the disc mate- 
rial. Moreover, since the orbit is described with respect to 
a varying rotating reference frame, all non-inertial terms in- 
volving the angular velocity, SI, and the angular acceleration, 
SI, of the coordinate system have to be taken into account. 
Restricting to those orbits coplanar with the disc midplane 
(8 = 7r/2 or z = 0), the equation of motion of the planet is 



. _ G(M.+M P ) 

Jp ~ |;cp| 3 p+ " 

-SI X (SI X x p ) -20xip-f!xx p . 



(7) 



We recall that, by working hypothesis, SI as well as SI are 
perpendicular to the disc midplane and produce a counter- 
clockwise rotation. The acceleration applied by the disc mat- 
ter to the planet is given by 



A P =G 



(x — x p ) dMo(x) 



1Mb (\x-Xp\2 + £ 2) 3 / 2 ' 
while the acceleration applied to the star is 



A*=G 



x dMu(x) 



M D 



(8) 



(9) 



In both cases the integration is carried out over the simulated 
disc mass, Mu (see section 2.4). 

Equation (8) contains the smoothing length in its denomi- 
nator. This expression for acceleration appears in the second 
term of equation (6). An acceleration with smoothing is ap- 
plied to the planet's motion in order to satisfy Newton's third 
law. 

2.3. Rotational elements of the reference frame 

The main aim of this paper is to study the exchange of 
angular momentum occurring between a migrating planet and 
disc material moving on U-turns of horse-shoe orbits. In order 
to accurately resolve the flow variables in this region by means 
of a local grid-refinement technique, the planet needs to move 
through the grid as slowly as possible. To achieve this, we 
worked in reference frames that rotate about the disc axis at 
a variable rate, SI = Sl(t). We then chose SI and SI so as to 
compensate for the fastest component of the planet motion, 
i.e., the azimuthal one. This is accomplished by calculating 
the total orbital angular momentum of the planet per unit 
mass, Ha, and then requiring that 



Ha = Xp x (Six Xp), 
Ha = XpX (A P - A,). 



(10) 
(11) 



These equations are to be solved with the additional require- 
ments that both SI and SI must be perpendicular to the plane 
of the orbit and produce a positive (i.e., counter-clockwise) 
rotation. Equations 10 and 11 constrain the angular veloc- 
ity and acceleration of the rotating coordinate system so that 
the planet trajectory reduces to a purely radial motion. In 
other terms, all of the planet's orbital angular momentum is 
conveyed to the rotation of the non-inertial reference frame. 
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If the orbital eccentricity remains close to zero during the 
system evolution, as we found in our simulations, then the 
planet's radial motion is only due to the disc gravitational 
torques. We denote the planet's semi-major axis as a = |x p | 
and the time-scale of this drifting motion as tm = a/\a\. The 
quantity N r Ar (or Nr AR) is the radial extent of the high- 
est refinement region and the time spent within this region 
is N r Ar/\a\ = N r (Ar/a) tm, which is on the order of 0.1 tm 
for the parameters used in the calculations. Numerical sim- 
ulations (e.g., Lubow et al. 1999; Nelson et al. 2000; Kley 
et al. 2001; D'Angelo et al. 2002) as well as analytical the- 
ories (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; 
Ward 1997) on disc torques suggest time-scales, tm, on the 
order of 10 4 periods. Therefore, with this method one can 
expect to track the planet and the coorbital regions, with the 
necessary numerical resolution, for hundreds of orbits. 

2.4. Physical parameters 

We performed two kinds of simulations: the first kind is 
dedicated to planets interacting with a low-mass disc and the 
second is dedicated to planets orbiting in a high-mass disc. 
In all of the calculations, the mass of the star, M* , represents 
the unit of mass whereas the initial semi-major axis of the 
planet's orbit, ao, gives the length unit. The unit of time is 
such that 1/to = \J G (Ml + M p )/a.Q. However, when it is 
necessary to convert quantities into physical units, we used 
hU = 1 M @ and a a = 5.2 AU. 

2.4-1- Parameters for low-mass discs 

In these models the simulated disc domain extends radially 
from 0.4 to 4.0 length units around the star and, azimuthally 
in angle, from to 27r. These simulations describe a disc 
of mass Md = 7.5 x 10~ 3 M* within the radial limits of the 
computational domain, which is equivalent to 0.01 Mq within 
26 AU of a 1 Mq star. In the case of 3D models, we simu- 
lated only the upper half of the disc between 80° < 9 < 90° 
and assumed mirror symmetry with respect to the midplane. 
The aspect ratio of the disc was fixed to H/r = 0.05. The 
overall initial surface density scales as r _1//2 and is axisym- 
metric. This would give an unperturbed disc surface density 
at the location of the planet of 76 gem -2 , but we included 
an initial gap along the planetary orbit that accounts for an 
approximate balance of viscous and tidal torques. One model 
was also run without an initial gap, in order to determine its 
influence on the results. In 3D models, the initial latitude 
dependence of the mass density is taken to be a Gaussian. 

We employed a constant kinematic viscosity, v, to account 
for the effects due to turbulence in the disc. In the units 
introduced above, we set v = 10 -5 that is also equivalent to 
Shakura & Sunyaev parameter a = 4 x 10~ 3 at the initial 
location of the planet. This choice is compatible with what 
was recently found in studies of embedded Jupiter-size bodies 
in discs with MHD turbulence (Papaloizou & Nelson 2003; 
Winters et al. 2003). However, we do not include the spatial 
variations in a consistent with the MHD results, nor the time 
fluctuations due to MHD turbulence (Nelson & Papaloizou 
2004). 

The planet mass is such that M p /M, = 10 -3 (i.e., one 
Jupiter-mass, Mj, for a one-solar-mass star). The planet 



starts on a circular orbit of semi-major axis ao = 1, which 
is kept static for a certain number of periods to allow the re- 
laxation of the system. This was done by setting to zero the 
terms (8) and (9), in equation (7), and activating them at the 
"release" time, t — t r i s . We used t r i s equal to either 100 or 
300 orbits. The migration rates were found to be insensitive 
to the release time (less than 10 per cent differences in rates), 
provided it is greater than 100 orbits. The azimuthal position 
of the planet remains constant throughout the computations 
(see section 2.3) and it is equal to (f> — (f> p = n. 

The smoothing length of the planet potential, e, in equa- 
tion (6) was chosen to be a fraction of the planet's Hill radius, 
Rh = a [A/p/(3A/,)] 1/3 = 0.069 a. We employed three differ- 
ent values: e = 0.4, 0.2, and 0.1 Rh, in order to study the 
effects of smoothing on the results. 

2-4-2. Parameters for high-mass discs 

When simulating planets embedded in a high-mass disc, 
we used parameters as similar as possible to those adopted 
by MP03, in order to have a direct comparison with their 
results. Therefore, in contrast to the previous settings, the 
radial extent of the disc and its aspect ratio were reduced to 
[0.4, 2.5] length units and 0.03, respectively. The simulations 
describe a disc of mass Md = 2.37 x 10 -2 M* within the radial 
limits of the computational domain, which is equivalent to 
24 Mj within 13 AU of a 1 M @ star. As in MP03, the initial 
surface density scales as r~ 3//2 and there is no initial gap. 
This gives an unperturbed disc surface density at the location 
of the planet of 653 gem -2 . The planet-to-star mass ratio 
is Mp/M* = 3 x 10~ 4 , roughly corresponding to a Saturn- 
mass object for M, = IMq. We again employed a constant 
kinematic viscosity v = 10~ 5 in dimensionless units. The 
planet was held on a static orbit (ao = 1) and released at t — 
t T \ B . For most of the 2D calculations the planet was released 
after 477 orbits, as done by MP03. For comparisons between 
2D and 3D models we could not afford the time required to 
run 3D calculations to 477 orbits, so we released the planet 
at 200 orbits. For a convergence test with high-resolution 
2D calculations we set t r is = 277 orbits. The value of the 
smoothing length was set to e — 0.3878 Rh (Rh — 0.046a), 
which is equal to 60 per cent of the local disc scale-height. 

3. Description of the numerical method 

The hydrodynamical equations (1) through (5) that de- 
scribe the evolution of the disc are solved numerically by 
means of a finite-difference scheme with directional operator 
splitting. The method is second-order accurate in space and 
first-order in time (Ziegler & Yorke 1997). The numerical res- 
olution of the regions around the planet is greatly enhanced 
by utilising a nested-grid technique (for details, see D'Angelo 
et al. 2002, 2003). Each subgrid level increases the resolution, 
with respect to the hosting grid, by a factor 2 in each direc- 
tion. Thus, the total gain in resolution for each added subgrid 
is 2 2 or 2 3 in 2D or 3D simulations, respectively. Subgrids are 
fully nested, i.e., each occupies a region of space completely 
contained inside the hosting grid. This implies that the num- 
ber of zones of any subgrid, along any direction, can be at 
most twice the number of zones of the hosting grid along that 
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Tabic 1: Grid system employed in low-mass disc models. 



Grid 


2D3G 


2D4G 


2D5G 


3D3G 


level 




N r xN^ 


N r x 


N R xN e x 


1 


243 x 455 


243 x 455 


243 x 455 


243 x 17 x 455 


2 


114 x 84 


114 x 84 


114 x 104 


114 x 24 x 84 


3 


114 x 84 


114 x 84 


134 x 124 


114 x 24 x 84 


4 




134 x 84 


154 x 144 




5 






174 x 164 





The linear resolution around the planet, averaged over all di- 
rections, on the level 1 is 1.45 x 10 -2 . This value decreases 
by a factor of 2^~ 1 - ) on a given level I. Thus, the grid sys- 
tems 2D3G and 3D3G resolve the flow around the Hill sphere 
of a Jupiter-mass planet with 19 grid zones per Hill radius, 
while the grid 2D5G achieves a resolution of 76 zones per Hill 
radius. 

direction. A point in space is handled by the highest resolu- 
tion grid (highest grid level) that covers that point. 

In order to test the behaviour of the nested-grid code for 
planetary migration calculations, we compared outcomes of 
models executed in a single-grid mode with those of the same 
models executed in a nested-grid mode with equal numeri- 
cal resolution. We always found an excellent agreement with 
discrepancies averaging ~ 10~ 3 per cent. Some of these com- 
parisons are reported in the Appendix. 

The equations of motion of the planet are solved in Carte- 
sian coordinates with a high-accuracy and fast hybrid algo- 
rithm. This involves a Bulirsch-Stoer method with an adap- 
tive time-step control (Press et al. 1992) and a standard 
fourth-order Runge-Kutta solver. Each hydrodynamics time- 
step At (constrained by the Courant-Friedrichs-Lewy stabil- 
ity criterion) is divided into substeps whose duration is dic- 
tated by the requirement that the local truncation error is 
always smaller than the chosen accuracy (10 -7 in these cal- 
culations). The maximum number of substeps allowed is set 
to 5000. If the integration time has not reached the value At 
after this iteration cycle, the remainder of the time-interval is 
integrated via a fourth-order Runge-Kutta method. Although 
this is a necessary precaution, the overall procedure actually 
requires only a few time-substeps of the Bulirsch-Stoer algo- 
rithm to complete the whole hydrodynamics time-step At, 
since the vector equation (7) has no singular points inside 
this integration interval. The orbit integrator was tested, 
over long-term evolutions, against both circular and eccen- 
tric Keplerian orbits. For a variety of values of £7 and f2, no 
deviations from the analytic solutions were found down to the 
machine precision. 

Disc gravitational forces given by equations (8) and (9) are 
considered to be constant over the whole time span At and 
are computed by summation of discretised quantities over the 
whole grid, always using densities from the subgrid with the 
highest resolution available. 

3.1. Numerical setup 

In a disc-planet interaction calculation, the largest spatial 
gradients of the flow variables develop around and inside the 



planet's Hill sphere. Over a distance of two Hill radii, the 
density can change by three or more orders of magnitude 
(D'Angelo et al. 2003; Bate et al. 2003) and the velocity field 
describes highly complex patterns (e.g., Lubow et al. 1999; 
Tanigawa & Watanabe 2002). In order to ascertain to what 
extent our numerical experiments depend upon the numerical 
resolution in this region, we performed a convergence study 
in all cases. To this aim, we set up a number of grid systems 
whose resolution in the coorbital regions ranges from 19 (or 
13, when M p /M« = 3 x 10~ 4 ) to nearly 76 (or 104) grid zones 
per Hill radius (see details in Tables 1 and 2). 

Most calculations were performed without allowing the 
planet to accrete. When accretion is permitted, mass is re- 
moved from around the planet within a tenth of its Hill radius. 
The removal of mass occurs on a time-scale on the order of a 
tenth of the orbital period. 

Boundary conditions at the inner radial border allow flow 
towards the central star, as naturally happens in viscous ac- 
cretion discs. The outer radial border is closed so that no 
inflow or outflow of material is permitted. At both radial 
edges of the disc, the flow is assumed to be Keplerian around 
the central star. This circumstance may occasionally lead to 
spurious, small-amplitude wave excitation at the outer edge 
of the disc, since material there has the tendency to orbit 
about the centre-of-mass of the system rather than around 
the central star (see discussion in Nelson et al. 2000). How- 
ever, the effects of such waves are not relevant since they do 
not propagate for a significant distance from the disc edge. 
In 3D models, reflective and symmetric boundary conditions 
are applied at the highest latitudes and at the midplane, re- 
spectively. The velocity field in the disc is initialised with a 
Keplerian circulation, corrected for the grid rotation. 

It was pointed out by Nelson & Benz (2003a) that when the 
planet moves through the grid cells, the smoothing length in 
equation (6) needs to be larger than half the linear dimension 
of the grid zone. This is required in order to avoid unphysical 
effects on the planet's trajectory due to close encounters with 
grid centres. In these calculations we used initial values of 
e that are at least 1.8 times the average linear size of the 
grid zone from which the planet is released. Furthermore, 
the series of convergence tests that we performed indicate 
that the ratios between e and the average grid zone size are 
large enough not to affect the outcome of the simulations (see 
section 4). 

4. Results of low-mass disc models 

Previous studies of migrating Jupiter- mass planets showed 
that the evolution of the orbital semi-major axis, a = a(t), 
is dependent upon the resolution with which hydrodynamics 
variables are discretised on the computing mesh (Nelson et al. 
2000; Nelson & Benz 2003a). As emphasised by Nelson & 
Benz (2003b), this issue becomes even more important when 
torques in the coorbital region are resolved, i.e., when the 
planetary gravitational smoothing lengths are a small fraction 
of i?H (see eq. [6]). The dependence of gravitational torques 
on the numerical resolution is crucial to assess the reliability 
of the outcomes. Therefore, we tackled this problem with 
a number of dedicated simulations, before investigating any 
possible physical effects of coorbital torques on the migration 
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Tabic 2: Grid system utilised in high-mass disc models. 



Grid 


2DlGb 


2D3Gb 


2D4Gb 


2D5Gb 


2D6Gb 


3D3Gb 


level 


N r X 


N r x N+ 


N r x 


X iV 




N R xN e x 7V 


1 


147 x 455 


147 x 455 


147 x 455 


147 x 455 


147 x 455 


147 x 17 x 455 


2 




114 x 84 


114 x 84 


114 x 84 


134 x 104 


84 x 24 x 84 


3 




114 x 84 


114 x 84 


114 x 84 


134 x 104 


84 x 24 x 84 


4 






114 x 84 


134 x 84 


164 x 104 




5 








164 x 104 


194 x 124 




6 










324 x 204 





The average linear resolution around the planet on the level I is 1.43 x 10 _2 /2" -1 '. Hence, the 
grid systems 2D3Gb and 3D3Gb resolve the flow around the Roche lobe of a 0.3 Afj planet with 
13 grid zones per Hill radius while the grid 2D6Gb achieves a resolution of more than 100 zones 
per Hill radius. We used the single-level grid system, 2DlGb, only for purposes of comparison 
with the calculations reported in MP03, 



of giant planets. 

4.1. A convergence study 

Convergence tests were carried out on each of the low-mass 
disc models described in Section 2.4.1. The resolution was 
progressively increased by employing the grid systems 2D3G, 
2D4G, and 2D5G (see Table 1). The last two grid systems, 
compared to the first, provide a linear resolution gain of a 
factor 2 and 4, respectively. For comparison purposes, we set 



3.9E-QS 7.9E-D5 1 .5E-D4 2.3E-04 3.1E-04 




Fig. 1. — Global surface density around a M p — 1 Mj planet 
orbiting in a low-mass disk (see section 2.4.1). The density is 
shown after 370 orbits, when the planet has migrated for 70 
orbits. In the linear grey-scale, at 5.2 AU, E = 10 -4 corre- 
sponds to 32.9 gem' 2 . 



the release time to t T \a = 100 orbits, except for the simulations 
concerning the accreting model that have t T i s = 300 orbits. 
The overall surface density from one of such calculations is 
displayed in Figure 1. 

Figure 2 shows that we achieved numerical convergence 
in all cases, with either accreting or non-accreting planets 
and with different values of e. In some panels, the out- 
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Fig. 3. — Surface density profile along the azimuthal position 
(4> = 4> P ) of a M p = lMj planet, after 100 orbital periods, 
orbiting in a low-mass disc. Different line types indicate mod- 
els with different values of the gravitational potential soften- 
ing: e — 0.1 Rh (solid line); e — 0.2 Rh (short-dash line); 
e = 0.4 Rh (dash-dot line); accreting planet (long-dash line). 
If a = 5.2 AU and M* = 1M Q , E = 10" 4 corresponds to 
32.9 gem" 2 . 
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Fig. 2. — Convergence tests regarding different configurations of Jupiter- mass models orbiting in a two-dimensional low-mass 
disc. Upper-left. Non-accreting planet with gravitational potential softening e = 0.4 Rh- Upper-right. Non-accreting planet with 
e = 0.2 Rh- Lower-left. Non-accreting planet with e = 0.1 Rh. Lower-right. Accreting planet with e = 0.1 Rh- The release time 
is equal to 100 orbits, except for the accreting model for which t T \ a = 300 orbits. Each panel shows how the semi-major axis, a, 
evolves when all torques from within the planet's Hill sphere are taken into account (upper curves) and when the contribution 
of those torques arising inside of 0.5 Rh from the planet are neglected (lower curves). See text for further details. 



comes produced by the two grid systems can be hardly distin- 
guished. The main numerical difficulty with the evaluation 
of gravitational torques arising from the coorbital region is 
related to the presence of large density gradients (Nelson & 
Benz 2003b). Moreover, the shorter the smoothing length, 
the larger such gradients are. Figure 3 shows that there is 
an order-of-magnitude difference between the density peaks 
of the models with e = 0.1 Rh and e = 0.4 Rh. For the 
e = 0.1 Rh case, using the grid system 2D3G would mean 



that e was resolved by less than 2 grid zones and, thus, the 
density gradients could be resolved too poorly. Therefore, we 
used the grid systems 2D4G and 2D5G for this case. 

We also investigated whether there are any differences be- 
tween simulations starting with or without an initial density 
gap (see section 2.4). Provided that the system is allowed to 
evolve for a sufficiently long time in order that the gap be- 
comes deep enough (ss 500 orbits), the migration behaviour 
is very similar to that of models initiated with a density gap. 
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Two sets of lines are displayed in each panel of Figure 2. 
These are intended to address the lingering question of the 
importance of torques exerted by matter residing deep inside 
the planet's Hill sphere (D'Angelo et al. 2003; Bate et al. 
2003). Thus, two configurations were simulated, differing only 
in whether or not torques within a radius (3Rh from the planet 
are included in the calculation of the gravitational force in 
equations (8) and (9). In one configuration, all torques are 
taken into account (i.e., j3 = 0). In the second configuration, 
the simulations were repeated neglecting the contribution of 
material lying inside the inner half of the Hill sphere (i.e., f3 = 
0.5). The choice (3 = 0.5 was made to avoid the region where 
the density gradient is largest (see Fig. 3) and which mostly 
contains material orbiting the planet before it is released (see 
Fig. 4). 

The streamlines in Figure 4 are constructed by integrat- 
ing the velocity field (u r ,u^) at an instant in time. Strictly 
speaking, this procedure is in error for a migrating planet 
(right panels), as it moves during the interval of integration. 
But provided the planet moves only a small fraction of its Hill 
sphere over the integration time, the streamlines obtained 
are reasonably accurate. This condition is satisfied for the 
streamlines plotted in this Figure. It is of course incorrect 
to ignore torques from within the Hill sphere since material 
can move into or out of it (see Fig. 3) and the angular mo- 
mentum associated with this mass flux is lost instead of being 
transferred to the planet's orbit. Therefore, migration rates 
obtained from configurations with /3 > are not fully consis- 
tent from a physical standpoint, unless one can assure that 
the neglected material is constantly and not temporarily or- 
biting the planet. 

In all cases we considered, the fi = calculations resulted 
in slower migration and, thus, these results appear as the 
upper curves in each panel of Figure 2. We note that even 
when torques coming from material within the Hill sphere 
are included, numerical convergence is still achieved. This 
last point turns out to be crucial when high-mass discs are 
considered (see section 5). 

4.2. Three-dimensional simulations 

The vertical stratification of the flow variables in (verti- 
cally isothermal) discs does not play an important role in de- 
termining the strength of Lindblad torques acting on Jupiter- 
mass planets, provided that Rh > H (Kley et al. 2001). How- 
ever, since the flow structure around the Hill sphere of the 
planet is fully three-dimensional (D'Angelo et al. 2003; Bate 
et al. 2003), the amount of angular momentum delivered by 
material in the vicinity of the planet may be affected by the 
vertical motion of the fluid. We attempted to investigate this 
issue by means of 3D calculations (grid system 3D3G), whose 
results are shown in Figure 5. As for 2D computations, the 
effects of torques exerted by material within the Hill sphere 
were measured by running models with (3 — and /3 = 0.5. 
In each panel of Figure 5, the evolution of the semi-major 
axis (solid lines) is compared to that obtained from 2D mod- 
els (dashed lines) having an analogous grid system (2D3G). 
We were unable to test for convergence of the 3D calculations 
due to computational limitations (simulations with a factor 2 
increase in linear resolution would have required around 5000 



CPU hours each). But since the 2D calculations were con- 
verged, we speculate that at the same resolution 3D calcula- 
tions are also converged because the density structure around 
and inside the Hill sphere is smoother in three dimensions. 

Figure 5 shows that the two- and three-dimensional results 
are similar in the non-accreting cases. The Hill spheres of 
non-accreting planets contain more material in 3D than they 
do in 2D (see Fig. 6). Therefore, it is reasonable to expect 
that the migration is slightly faster in three dimensions. 

The situation appears more complex in the accreting case, 
for which migration is slower in 3D than in 2D if /3 = 0.5, 
but it is faster if /3 = (Fig. 5, right panel). Since the mass 
inside the Hill sphere is nearly the same in the two geome- 
tries (Fig. 6, bottom panel), we ascribed this discrepancy to 
the strong spiral waves that occur in the two-dimensional ac- 
creting calculations (Lubow et al. 1999; D'Angelo et al. 2002) 
which are much weaker in three-dimensions due to the possi- 
bility of vertical motions. Strong spiral waves do not develop 
when e is a fair fraction of _Rh- The non-accreting models 
with e = 0.4 and 0.2 Rh present a nearly featureless density 
structure close to the planet in both geometries, hence the 
similarity of the migration rates. 

Finally, we note that non-accreting 3D migration rates 
with e = 0.4 or 0.2 _Rh are very similar. This suggests that no 
large variations should be expected if smaller values of e were 
employed, provided that a sufficiently refined mesh is utilised 
(A_R = RAtf) <C e). The same conclusion seems to be valid in 
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Fig. 6. — Mass enclosed within the Hill sphere of a M p = 
lMj planet, as measured in 2D (circles) and 3D (diamonds) 
calculations. Top. Non-accreting model with e — 0.4 _Rh- 
Centre. Non-accreting model with e = 0.2 Rh. Bottom. Ac- 
creting model with s = 0.1 Rh. Whether or not a gap is im- 
posed on the initial density structure, the amount of material 
within the Hill sphere tends toward the same value. 
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Fig. 4. — Surface density and streamlines around Jupiter-mass, non-accreting planets orbiting in a low-mass disc (Md = 
0.01 Mq within 26 AU). The results were obtained from computations executed with the grid system 2D4G (linear resolution 
Rh/ = 38) and softening e = 0.4_Rh (top) and 0.1 Rh (bottom). The panels illustrate the situation at the release time t — 300 
orbital periods (left) and 50 orbits later (right), while the planet is migrating. The grey-scale is logarithmic and, at 5.2 AU, 
S = 10~ 3 corresponds to 329 gem -2 . The two streamlines closest to the planet start from distances of w 0.36 and ~ 0.5 Rh, 
respectively. 



two dimensions, as discussed more quantitatively in the next 
section. 

4.3. Migration rates: a quantitative analysis 

So long as the semi-major axis does not change signifi- 
cantly (i.e., it remains of the same order of magnitude), the 
migration of a Jupiter-mass planet roughly follows an expo- 
nential decay (Nelson et al. 2000). We assume that even when 



the action of torques arising from corotation regions and from 
material orbiting the planet are included, the evolution of a 
can also be described by an exponential decay law 

a(t) = a e- (t - trls)/TM , (12) 

for times t > t T i s . We take the migration time-scale tm = 
a/\a\ to be a constant over the simulated time-interval of the 
actual planet's migration (between 40 and 100 orbits). This 
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Fig. 5. — Three-dimensional simulations of Jupiter-mass models (grid system 3D3G) orbiting in a low-mass disc, compared 
to analogous two-dimensional models (grid system 2D3G). Left. Non-accreting planet with gravitational potential softening 
s = 0.4 _Rh- Centre. Non-accreting planet with e = 0.2 Rh- Right. Accreting planet with e = 0.1 Rh- The release time is equal 
to 100 orbits for the non-accreting models and 300 orbits for the accreting models. In each panel, the slower migration occurs 
with the configuration executed with /3 — 0, whereas the faster migration occurs with the configuration run with /3 = 0.5. 



simple parameterisation of a = a(t) is very useful because 
tm can be directly connected to the acting torques. In fact, 
if the orbit eccentricity is negligible then the conservation of 
the orbital angular momentum leads to the relation 



2T ■ n 



(13) 



in which the vector T denotes the total external torque. This 
expression is commonly used to evaluate a from the vertical 
component of T (which we simply indicate as T), when a 
planet moves on a static orbit (i.e., it is not allowed to mi- 
grate) . 

We obtained estimates of tm for all of the 2D models (grid 
system 2D4G) by performing a linear least-mean-squared fit 
of the relation ln(a/eto) = — (t — tris)Ato- The results are 
listed in the second and third columns of Table 3, labelled 
as "moving" migration time-scale. The relative error on each 
estimate is at most 10 -3 and only for this reason tm is given 
with three significant digits. Discrepancies between estimates 
computed from 2D and 3D non-accreting models are below 
~ 10 per cent. Since the accreting models present a more 
significant discrepancy, tm is also reported for the simulations 
in three dimensions. 

Including the effect of matter orbiting the planet tends to 
slow down its inward drifting motion, regardless of the em- 
ployed disc geometry, as clearly indicated in Figure 5. The 
comparison between the (3 = and = 0.5 migration time- 
scales shows that the torque from this material can be compa- 
rable to that from corotation and Lindblad resonances. The 
total (positive) torque produced inside the inner half of the 
Hill sphere is 



Tks = T — T~lc oc 



1 



r M [/3 = 0] t m [/3 = 0.5] : 



(14) 



whereas the magnitude of the total (negative) torque ex- 
erted from the rest of the disc (i.e., Lindblad and corotation 



torques), |7lc|j is proportional to 1/tm[/3 = 0.5]. Hence, the 
ratio between the two contributions is 



Ths 
l^c| 



tm[/? = 0.5] 
tm\P = 0] ' 



(15) 



Entries in the second and third columns of Table 3 indicate 
that in the model with softening e = 0.4 Rh the material close 
to the planet accounts for a relatively small contribution (22 
per cent). However, shorter smoothing lengths dramatically 
increase the torque ratio, which becomes greater than 60 per 
cent in the non-accreting models with s = 0.2 Rh and 0.1 Rh. 
A similar ratio between torques is obtained in the 3D accret- 
ing model. 

These migration time-scales can be compared with the 
Type I (no gap, resonant) time-scale of about 4 x 10 2 or- 
bits in 2D and 6 x 10 2 orbits in 3D (Tanaka et al. 2002) and 
the Type II (viscous) time-scale 2a 2 /(3^) ~ 10 4 orbits. 

Note that the 2D migration rates tend to converge as e is 
decreased. In particular, the migration rates for e = 0.2 Rh 
and 0.1 Rh differ by less than 10 per cent. 

4.4. Comparison of migration rates of static and 
migrating planets: the Jupiter-mass case 

We examined whether the torque exerted on the planet 
by the disc material is influenced by the radial motion of the 
planet. As discussed in Section 1, the motion of the planet 
might be able to affect the coorbital torques and therefore 
the migration rate. In order to test this hypothesis, we com- 
puted the total torque acting on the planet during the last 
ten orbital periods before it was released. This was done for 
both j3 = and (3 — 0.5 configurations. Since no angular 
momentum is actually extracted from or added to the plan- 
etary orbit, which thus remains static, we shall refer to such 
torques as static torques. The migration time-scales listed in 
the two right-most columns of Table 3 were obtained from the 
average static torques by using equation (13). In Table 3 we 
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compared these "static" migration time-scales, t"m, with the 
migration time-scales, tm, measured from the moving planet 
calculations. In all cases, there is close agreement between the 
static and moving migration time-scales. These results show 
that under these circumstances of disc and planetary masses, 
there is no strong dependence of the torques on whether plan- 
ets are on fixed orbits or allowed to migrate. 

5. Results of high-mass disc models 

The Type II migration rate depends only on the viscous 
time-scale of the disc near the location of the planet and is in- 
dependent of the disc density, provided that the gap is devoid 
of material. Yet gaps are generally not completely cleared and 
the Type II time-scale prediction does not take into consider- 
ation the angular momentum exchanged between the planet 
and the "gap" material. Some of this material travels on 
horse-shoe orbits, while other material circulates within the 
planet's Hill sphere. The angular momentum delivered in ei- 
ther case may play a major role in planetary migration (see, 
e.g., Masset 2001) and it is proportional to the local mass 
density. In fact, MP03 recently claimed that there exists a 
critical mass (when the material around the planet is more 
massive than the planet), beyond which a runaway migration 
process sets in. 

We ran simulations of Saturn-like bodies (M p = 0.3 Mj) 
embedded in a disc as massive as 24 Mj inside 13 AU. The 
annular region within 2 _Rh from the planet is initially 7.5 as 
massive as the planet. Nonetheless, the aspect ratio is small 
enough (H/r — 0.03) so the thermal condition for gap for- 
mation, Mp/M* > 3 (H/r) 3 (e.g., Lin & Papaloizou 1993), 
is fulfilled and therefore the migration might be within the 
Type II regime, although the gap is not completely cleared 
as can be seen in Figure 7. In these cases, with massive discs 



Table 3: Comparison of static and moving migration 
time-scales for a Jupiter-mass planet in a low-mass disc. 



Moving 



Static 



e 


(3 = 


(3 = 0.5 


/3 = 


/3 = 0.5 


0.4 i? H 


9.94 x 10 3 


7.75 x 10 3 


1.0 x 10 4 


7.8 x 10 3 


0.2 i? H 


1.76 x 10 4 


6.34 x 10 3 


1.8 x 10 4 


6.2 x 10 3 


0.1 i? H 


1.51 x 10 4 


5.71 x 10 3 


1.8 x 10 4 


5.7 x 10 3 


0.1 


5.97 x 10 4 


4.86 x 10 3 


4.8 x 10 4 


4.2 x 10 3 


0.1 Rh 1 


1.54 x 10 4 


5.77 x 10 3 


1.5 x 10 4 


5.4 x 10 3 



t 2D accreting model. * 3D accreting model. 

The migration time-scales labelled as "moving" refer to the 
time-scale, tm, in equation (12) and were computed as ex- 
plained in Section 4.3. They are expressed in units of initial 
orbital periods, i.e., 11.9 years if ao = 5.2 AU. One-standard 
deviation uncertainties for these estimates range from 1 to 10 
orbits. See Section 4.1 for an explanation of configurations 
/3 = and — 0.5. Migration time-scales labelled as "static" 
were determined from equation (13) by employing torques av- 
eraged over the last ten orbits before the release time. Com- 
putations were executed with the grid system 2D4G. 



and small aspect ratios, very large density gradients develop 
inside the Hill sphere. Therefore, it is especially important to 
investigate the dependence of the results on numerical resolu- 
tion. We achieved convergence for the flow outside of the Hill 
sphere by using numerical resolutions of order 13 grid zones 
per Hill radius. However, in order to accurately determine 
the contributions to the migration rate from material inside 
the Hill sphere, resolutions higher than 52 grid zones per Hill 
radius are necessary. 

5.1. Convergence tests 

As mentioned in Section 2.4.2, the model setup and the 
disc parameters were chosen to match as closely as possible 
those in MP03. The smoothing length was 60 per cent of the 
local disc thickness, H, (i.e., e = 0.3878 Ru), the planet was 
non-accreting, and t T \ s = 477 orbits. We performed a calcu- 
lation using a single-level grid (2DlGb, see Table 2) aimed 
at reproducing the resolution used by MP03 (Ar/Ru — 0.3 
and Ar/e ~ 0.8). We then performed a convergence test us- 
ing different numerical resolutions, as provided by the grid 
systems 2D3Gb, 2D4Gb, and 2D5Gb (see Table 2). An addi- 
tional convergence test, involving the grid system 2D6Gb, is 
discussed in Section 5.1.1. 

The left panel of Figure 8 shows the outcomes of the tests 
for the evolution of the semi-major axis concerning the con- 
figuration with (3 = 0. The dot-dash line in this panel rep- 
resents the result from the single-grid computation 2DlGb, 



-4.44 -3.H2 -3.21 -2.BQ -1.99 



^ D - 



-1 - 




Fig. 7. — Global surface density around a M p = 0.3 Mj (non- 
accreting) planet orbiting in a 0.024 Mq disk. The density is 
displayed at t = 550 orbits, when the planet has migrated for 
about 70 orbits. The grey-scale is logarithmic and 10~ 3 cor- 
responds to 329 gem" 2 , at 5.2 AU. The average gap density 
is 2.5 x 10~ 4 or 82 gem" 2 . 
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Fig. 8. — Computations of a Saturn-mass planet orbiting in a high-mass disc: convergence tests. The gravitational potential 
softening, e, is 60 per cent of the local disc scale-height. The release time is equal to 477 orbits. The left panel shows the 
evolution of a when all torques are taken into account (i.e., (3 = 0). The right panel shows how the evolution of the semi-major 
axis proceeds when the contribution of those torques arising from inside the Hill sphere (i.e., /3 = 1.0) are neglected . The 
dot-dash line (labelled as 2DlGb) refers to a calculation executed with the same numerical resolution as in MP03. 



which should be compared to the model labelled as Ss in Fig- 
ure 2 of MP03. Given the remarkable agreement between our 
and their outcome, we are confident that we reproduced the 
same physical and numerical conditions for runaway migra- 
tion. Yet, computations repeated with finer and finer res- 
olutions gave smaller and smaller migration rates which, as 
displayed in Figure 8 (left panel), failed to converge. The gain 
in linear resolution achieved (over the single-grid simulation) 
with the employed grid systems ranges from 4 (2D3Gb) to 
16 (2D5Gb). In the highest resolution models, there are 52 
grid zones per Hill radius. Comparing the short-dash and 
dot-dash curves in left panel of Figure 8, one realises that 
the average migration speed obtained over the first 25 or- 
bits with the grid system 2D3Gb is only half (in physical 
units, (d) ~ -5 x 10 -3 AUyr" 1 ) of that in MP03. Calcu- 
lations executed with the grid systems 2D4Gb and 2D5Gb 
give even lower migration speeds of (a) ~ —5 x 10~ 4 and 
— 1.4 x 10~ 4 AUyr -1 , respectively. 

While there is a factor of 10 decrease in disc torques acting 
on the planet in going from grid systems 2D3Gb to 2D4Gb, 
this factor reduces to 3.6 when the two most refined grid 
systems are considered. Yet, from the behaviour of semi- 
major axis evolution shown in the left panel of Figure 8, we 
cannot determine whether it is converging. To assess this 
point we employed the grid system 2D6Gb (see section 5.1.1) 
which indicates that the solid line in Figure 8 is basically a 
converged evolution. 

The right panel in Figure 8 shows the semi-major axis 
evolution from the same calculations as in the left panel but 
executed with (5 = 1.0, i.e., excluded torques arising inside the 
planet's Hill sphere. As before, this choice of /3 was made to 
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Fig. 9. — Mass enclosed within the Hill sphere of a M p = 
0.3 Mj planet, as measured from simulations with increasing 
resolutions: single-level grid 2DlGb (asterisks); grid system 
2D3Gb (diamonds); grid system 2D4Gb (squares); grid sys- 
tem 2D5Gb (circles). 
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exclude the region around the planet with largest density gra- 
dients as well as largest torque densities. Clearly, numerical 
convergence was readily achieved with this configuration. The 
migration time-scale, obtained from a least-mean-squared fit 
to the data (see section 4.3), is tm = 493 orbits. Furthermore, 
outcomes of simulations executed with (3 = 0.75 attained con- 
vergence at almost the same rate of migration as with (3 = 1.0. 
Therefore, we conclude that the material close to the planet 
must be held responsible for the non-convergence of the (3 — 
configuration in the left panel of Figure 8. That is, the torque 
arising from within the Hill sphere converges very slowly with 
increasing resolution. Despite the fact that the amount of ma- 
terial inside the planet's Hill sphere increases as the grid res- 
olution is raised (see Fig. 9), the resulting migration rates or 
net torques are actually smaller. Figure 10 shows the surface 
density near the planet at an advanced time, shortly before it 
is allowed to migrate. We note that the mass near the planet 
seems to be converging at the highest resolutions, but conver- 
gence is not yet formally achieved. This Figure also implies 
that most of the material is piled up very close to the planet. 
We measured that 80 per cent of the mass contained inside 
the Hill sphere is concentrated within a distance of 0.2 Ru 
from the planet. 

The mass build up within the Hill sphere appears to sug- 
gest that disc self-gravity may be dynamically important. 
However, this may not actually be the case. A simple cal- 
culation of a viscous disc that accretes at the typical rates of 
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Fig. 10. — Surface density profile at the azimuthal posi- 
tion of a planet with M p /M* = 3x 10~ 4 , after 450 orbital 
periods (the planet starts migrating at t T i s = 477 orbits), 
orbiting within a high-mass disc. Different line types re- 
fer to computations performed with different grid systems: 
2D5Gb {solid line); 2D4Gb (long-dash line); 2D3Gb (short- 
dash line). If ao = 5.2 AU and M, = 1 AT©, £ = 10~ 2 is equal 
to 3.29 x 10 3 gcm~ 2 . 



10 -8 Mq per year suggests that it is likely not self-gravitating 
(the value of the Toomre parameter Q is much greater than 
unity for the parameters in this paper). In the simulations 
presented here, the mass build up is concentrated in a re- 
gion of order the smoothing length (see Fig. 10). Within that 
radius, further inward viscous accretion is artificially slow be- 
cause the gravitational potential of the planet tends to enforce 
rigid rotation (see second term in eq. [6]). In addition, the 
boundary condition of no accretion on the planet prevents the 
accumulated gas from being removed from the simulation. As 
such, much of the gas accumulated within a smoothing length 
represents material that is incorporated by the planet, rather 
than residing in the disc. 

Although the configuration with (3 = 1.0 provides nu- 
merically converged migration rates, one has to be wary of 
their physical meaning. Figure 11 illustrates that before the 
planet is released (top-left panel) material on horse-shoe or- 
bits passes through the Roche lobe as close to the planet as 
~ 0.5 Rh- Yet, if the planet starts rapidly migrating this 
picture is bound to change. The top-right panel shows a 
snapshot after 20 orbits from the release time, as the planet 
radially moves at a speed a ~ —1.5 x 10 _4 AUyr _1 (con- 
figuration (3 — 1.0 and grid system 2D5Gb). The situation 
appears less symmetric than before the release and the flow 
structure within the Hill sphere has been altered by the rapid 
planetary motion. As a reference, we also show in the bottom 
panel what happens when all torques are consistently taken 
into account (/? = 0). We calculated the torques arising from 
the Hill sphere in the situation depicted in the top-right panel 
and we found that they are three times as large (and more 
positive) as those exerted, at the same time, in the configura- 
tion (3 = (bottom panel). This difference may indicate that 
the faster motion in the (3=1 case has artificially changed 
the density distribution inside the Hill sphere and thus the 
circulation in the coorbital region. 

The reason for the very fast migration rate, measured at 
the lowest resolution (single-level grid 2DlGb), can be un- 
derstood by examining the two dimensional linear map of the 
torque density magnitude in the left panel of Figure 12. The 
plot describes the situation after 19 orbits from the planet's 
release, when it is migrating inwards at an average rate of 
roughly 10 -2 AUyr -1 . The map clearly shows how the poor 
resolution (the grid zone size is indicated by the shaded pix- 
els) cannot properly handle the large torque gradients within 
the Hill sphere and produces a very large differential torque. 
This resolution effect led to the vastly different migration 
time-scales between the lowest and highest curves in the left 
panel of Figure 8. A cut of the torque density magnitude, 
through the planet's radial position, is shown in the right 
panel of Figure 12 for both the computation executed with 
the grid 2DlGb (solid line) and that executed with the high- 
resolution grid system 2D5Gb (dashed line). The dashed-line 
profile was rescaled so that the maximum values were similar 
to those of the solid-line profile. The filled circles represent 
the actual data. The low-resolution torque density is highly 
asymmetric. The two maxima alone exert a negative torque 
that would result in a migration time-scale of 80 orbits. The 
large mismatch between the torque density extrema is not ob- 
served in high-resolution model, in which their opposite sign 
contributions nearly cancel each other. 
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Fig. 11. — Surface density and streamlines around a 
M p = 0.3 Mj, non-accreting planet orbiting in a high- 
mass disc (M D = 0.024 M Q within 13 AU). The dashed 
line indicates the Roche lobe and the crosses mark the 
position of the LI and L2 Lagrange points. The top- 
left panel refers to the time t = 450 orbital periods 
(i.e., before it is released) whereas the top-right panel 
displays the situation at t = 497 orbits with the con- 
figuration p = 1.0 (its instantaneous radial speed is 
d ~ —1.5 X 10 -4 AUyr - 1 ). The bottom panel refers to 
the same time but when all torques are taken into ac- 
count (i.e., (3 = 0). The grey-scale is logarithmic and, at 
5.2 AU, 10~ 3 corresponds to 329 gem -2 . These results 
were obtained from computations executed with the grid 
system 2D5Gb (linear resolution Ar/Ru ~ 2 X 10 -2 ). 
As in Figure 4, the streamlines in the bottom panel 
do not account for the planet's motion because of the 
small d. A more strict procedure was instead employed 
to calculate the streamlines in the top-right panel. In 
this case we integrated the velocity field (u T — d, u^), 
where d is the instantaneous radial speed of the planet. 



The differences between the two highest resolution calcu- 
lations discussed here are less evident and require some dis- 
cussion. Two-dimensional logarithmic maps of the magni- 
tude of the torque density for such models are shown in the 
top panels of Figure 13. They were obtained from the com- 
putations with the grid systems 2D4Gb (left) and 2D5Gb 
(right). Both maps describe the situation 60 orbits after the 
planet's release. The torque density is positive on the side 
leading the planet, cj> > <j> p , and negative on the opposite 
side (<f) < (j> p ). As clearly indicated in the Figure, the torque 
density within the inner half of the Hill sphere is orders of 
magnitudes larger than it is anywhere else in the surrounding 
region and, therefore, in the whole disc. This is the reason 
why the migration speed is so susceptible to the torques ex- 
erted within the planet's Hill lobe. Any mismatch between 



the positive (</> > <f> p ) and negative (ci > (f> p ) contributions 
can produce a very large net (either positive or negative) 
torque acting on the planet. From the top panels in Fig- 
ure 13, the torque density magnitude appears rather symmet- 
ric with respect to the direction <j> = cj> p . This is clear from 
the bottom-left panel, where cuts through the planet's radial 
position are compared for the two grid systems. The solid 
line corresponds to the more resolved model. Nevertheless, 
the results shown in the bottom-left panel of Figure 8 imply 
that the torque exerted by the Hill sphere is more positive 
(i.e., greater than zero and larger) in the higher resolution 
model (grid system 2D5Gb) than it is in the lower resolution 
one (grid system 2D4Gb). Indeed, this effect is highlighted 
by the ratio between the two torque density cuts (solid to 
dashed profile, that is higher to lower resolution results) re- 
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Fig. 12. — Left. Absolute value of the torque density (linear scale) close to a planet with M p /M« = 3 x 1CP 4 and orbiting in a 
high-mass disc. The dashed line indicates the Roche lobe and the crosses mark the LI and L2 Lagrange points. This calculation 
was executed with the single-level grid 2DlGb (Ar ~ 0.3 -Rh). Shaded pixels represent the actual size of the grid zones. The 
torque distribution is illustrated at t — 496 orbital periods while the planet is migrating at an average rate (d) ~ — 1CP 2 AUyr -1 . 
The torque density is negative when cj> < (f> p — tt and positive when cj> > cj> p . It is evident that at such low resolution there is an 
unbalanced inward torque that is not observed in higher resolution calculations (see top panels of Fig. 13). Right. The solid line 
shows the profile of the absolute value of the torque density (shown in the left panel) through the radial position of the planet. 
The filled circles indicate the positions of the data. The dashed line refers to an analogous profile from the highest resolution 
calculation (grid system 2D5Gb), which was rescaled so that its maximum value was 1.3 x 10 -6 . 



ported in the bottom-right panel of Figure 13. The important 
thing to note is that the curve is asymmetric, with respect to 
the direction <j> = cj> p , towards the the outer parts of the Hill 
sphere \cj> — (j) p \ > 0.02 ~ OARu/a. This means that the mis- 
match between the positive (</> > <f> p ) and negative (<f> < 4> p ) 
torques arising from the region \<j> — (f> p \ > 0.02 produces a 
net positive torque that is greater in the higher resolution 
model than it is in the lower resolution model. Most of the 
asymmetry, and therefore the discrepancy between the two 
computations, must be confined to the region enclosed be- 
tween roughly 0.4 Rh and 0.75 Rn from the planet because 
convergence tests executed with the configuration j3 = 0.75 
gave the same migration behaviour for the two models. 

We also performed 3D simulations, using the grid system 
3D3Gb (see Table 2), yet no appreciable differences from 2D 
calculations executed with the grid system 2D3Gb were ob- 
served. This was easily predictable, given the large smoothing 
length adopted in these models and the very small the disc 
thickness. 

5.1.1. A converged migration rate 

In order to evaluate how close to convergence the orbital 
evolution given by the grid system 2D5Gb is (see Fig. 8, left 
panel), we made a final attempt and ran a model with the 
grid system 2D6Gb (see Table 2), which resolves the Hill ra- 
dius with about 104 grid zones. However, we could not run 
a complete model as those in Section 5.1. In fact, evolving a 



model for about 550 orbits with such a grid system would have 
required around 8000 CPU hours. We only had the compu- 
tational resources to run this particular model for about 292 
orbital periods. Therefore we set the release time to t T i s = 277 
orbits and let the planet migrate for about 15 orbits. 

To carry out a consistent comparison, we performed cal- 
culations with the grid systems 2D4Gb and 2D5Gb imposing 
the same release time. The results are shown in Figure 14. 
Despite the short time over which the planet actually mi- 
grated, the highest resolution model (solid lines) provided 
evolutions of both a (top panel) and d (bottom panel) that 
are in very good agreement with those computed with the grid 
system 2D5Gb (long-dash lines). This implies that the rates 
of migration obtained with the latter grid system (2D5Gb) 
can be considered as converged rates. This also indicates 
that in order to accurately compute torques from within and 
around the Hill sphere, in a configuration as described in Sec- 
tion 2.4.2, linear resolutions on the order of 52 grid zones per 
Hill radius are required. 

It has to be emphasised that, while for Jupiter-mass plan- 
ets in low-mass discs torques were converged with respect to 
both the numerical resolution and the smoothing length of 
the planet's potential (see section 4), in the present case we 
only examined convergence with respect to the numerical res- 
olution. 
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Fig. 13. — Top. Absolute value of the torque density (logarithmic scale) and contour lines around and inside the Hill sphere of 
a planet with M p /M, = 3 x 10~ 4 , orbiting in a high-mass disc. The snapshots are taken at t = 537 orbital periods, i.e., after 
the planet has migrated for 60 orbits. The left panel refers to the calculation executed with the grid system 2D4Gb while the 
right panel refers to that executed with the grid system 2D5Gb. The torque density is negative when <f> < <j> p and positive when 
<j> > 4> p . Note that the torque density at the planet position is not exactly zero because the planet does not sit on the centre of 
a mesh zone, where density torque is computed. Bottom. The left panel shows the cut of the torque density magnitude through 
the planet's radial position for the models shown in the top panels. The solid line represents the outcome of the higher resolution 
simulation (2D5Gb) whereas the dashed line corresponds to the lower resolution simulation (2D4Gb). The ratio between these 
two curves (higher to lower resolution calculation) is displayed in the right panel. The profile appears asymmetric (with respect 



to 



0) only towards the outer regions 



> 0.02 ~ 0.4i? H /o. 



5.2. Comparison of migration rates of static and 
migrating planets: the Saturn-mass case 

We performed the same type of comparison, as done above, 
for the torques acting on a static and migrating planet. We 



considered both the configurations f3 = and f3 — 1.0. These 
results were obtained from the calculation run with the grid 
system 2D4Gb (26 grid zones per Hill radius). Recall that 
with /3 = 1.0, complete numerical convergence was attained 
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with 13 grid zones per Hill radius and therefore the same re- 
sult is produced by the model executed with the grid system 
2D5Gb. With j3 = 0, convergence was presumably obtained 
only with this last grid system. Nonetheless, it is still of inter- 
est to investigate if migration times-scales depend on whether 
the planet is allowed to migrate or kept on a fixed orbit, with 
the grid system 2D4Gb, since it yields a larger migration rate. 

Using the torques measured during the last ten orbits be- 
fore the planet is released, by the same procedure outlined 
in Section 4.4, we obtained static migration time-scales of 
t m = 770 initial orbits for /3 = and t m = 540 initial or- 
bits for /3 — 1.0. Both values are remarkably close to those 
obtained when the planet is allowed to move: tm = 773 and 
tm = 493 initial orbits, respectively (see Fig. 8). 

6. Discussion and Conclusions 

We calculated the migration rates of planets embedded in 
discs. These calculations were performed in two and three 
dimensions, using a reference frame that corotates with the 
planet and a nested-grid code that can provide high resolution 
close to the planet while it migrates. The models span a vari- 
ety of smoothing parameters for the potential and a variety of 
grid resolutions. Both accreting and non-accreting boundary 
conditions near the planet were considered. We were espe- 
cially interested in whether torques from the coorbital region 
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Fig. 14. — Semi-major axis evolution (top) and migration 
speed (bottom) obtained from the grid system 2D6Gb (solid 
line) versus those calculated with grid systems 2D5Gb (long- 
dash line) and 2D4Gb (short-dash line). The units of a are 
initial Hill radii per orbit. In these simulations the release 
time was equal to t T \ s — 277 orbits (see text). All torques are 
taken into account (i.e., /3 = 0). This additional convergence 
test shows that the migration behaviour given by the grid 
system 2D5Gb is fundamentally a converged one. 



can lead to runaway migration, as reported in MP03. 

In the case of a Jupiter-mass planet embedded a low-mass 
disc (Md = 0.01 Mq within 26 AU), the planet opens a gap in 
which there is some flowing material. Numerical convergence 
was readily obtained (see Fig. 2). The torques arising from 
within the Hill sphere do make a contribution to the torque 
(up to 60 per cent), but always in the sense of reducing the 
migration rate. The migration time-scales are numerically of 
order the Type II migration time-scale 2 a 2 / (3 v) ~ 10 4 orbits 
(see Table 3) and much longer than the Type I time-scale of 
about 5 x 10 2 orbits. 

In the case of a Saturn-mass planet embedded in a high- 
mass disc (Md = 0.02 Mq within 13 AU), the planet opens 
a less clean gap and is much more susceptible to the larger 
amount of material that resides in the coorbital region. Nu- 
merical convergence in this case was much more difficult to 
achieve when torques from within the Hill sphere were in- 
cluded. Convergence was more easily obtained when consider- 
ing only torques from outside the Hill sphere (see Fig. 8). The 
reason is that the mass of gas flowing within the Hill sphere 
is larger than the mass of the planet. Any inaccuracies in the 
density structure near the planet (e.g., due to finite resolu- 
tion) can lead to strong net torques (see Fig. 12 and Fig. 13). 
Although numerically converged, migration rates that do not 
account for torques from within the Hill sphere (or a large 
fraction of it) are artificially large (compare curves for grid 
system 2D5Gb in the left and right panels of Fig. 8). This also 
affects the flow structure around the planet and in the coor- 
bital region (see Fig. 11). With increasing resolution, the gas 
mass within the Hill sphere increases, yet the migration rate 
decreases. In the case that the resolution was about equal to 
the smoothing length, which was 0.39 Rh, the migration rate 
was very high, comparable to the Type I migration rate (as 
also found by MP03). However at the highest resolution we 
applied with a release time of 477 orbits, which was sixteen 
times higher (in terms of linear resolution), the migration rate 
dropped dramatically by more than two orders-of-magnitude. 
Discrepancies between higher resolution calculations are more 
subtle and arise from the outer half of the Hill sphere (see 
Fig. 13). A calculation based on the grid system 2D6Gb, 
for which Rn/Ar = 104, indicates that accurately describing 
torques from around and inside the Hill sphere requires res- 
olutions of at least 52 grid zones per Hill radius. At highest 
resolution, the migration time-scale is about 3 x 10 3 orbits 
(see left panel of Fig. 8), somewhat shorter than the Type II 
time-scale. But the process is unlikely to be simply described 
by Type II migration. 

We calculated the torques exerted by the disc on plan- 
ets whose orbits are fixed and used these to obtain migra- 
tion time-scales. Comparing these time-scales to those ob- 
tained by releasing the planet and allowing it to migrate freely 
through the disc, we found no significant difference in the mi- 
gration time-scales. This argues that the corotation torques 
are not greatly affected by the radial drift of a planet. 

In summary, the migration rates for planets that open im- 
pure gaps (in which some material flows) are substantially 
smaller than Type I rates and do not seem to be simply de- 
scribed by Type II migration. Torques arising near the planet 
can be important, but do not appear to have a dramatic effect 
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in raising the rates. Resolution is key to obtaining accurate 
torques. 
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Numerical tests 



The purpose of this Appendix is to demonstrate the re- 
liability of the nested-grid technique when it is applied to 
disc-planet interaction calculations and, more specifically, to 
planetary migration. The capabilities of this technique in the 
context of astrophysical fluid-dynamics modelling have been 
addressed by a number of authors (e.g., Ruffert 1992; Yorke, 
Bodenheimer & Laughlin 1993; Yorke & Kaisig 1995; Ziegler 
& Yorke 1997, and references therein). We therefore concen- 
trate on specific test computations that closely concern the 
application we did in this paper. The tests that we report 
here were done in two dimensions only to avoid excessively 
long computing times. 

We set up a model of a Jupiter-mass planet (M p /M* = 
1CT 3 ) in a massive disc with Mr, = 1.1 x 1CT 2 M* (inside 
13 AU of a 1 Mq star) and whose aspect ratio is H/r = 0.03. 
The initial surface density drops as r -1 / 2 but it also includes 
a theoretical gap along the planetary orbit, as done for most 
of the Jupiter-mass models discussed so far. The adopted 
viscosity prescription was the same as that chosen for all the 
other calculations (see section 2.4.1). The radial extent of the 
computational domain ranges from 0.4 to 2.5 length units and 
reflective boundary conditions were applied to both edges of 
the disc to enforce mass conservation (the planet was not ac- 
creting). The planet's orbit was initialised with ao = 1 (and 
with zero eccentricity) and it was kept steady for the first 
orbital period (i.e., t T i s — 1). The reference frame was set to 
rotate at a variable rate = £l(t) so that cj> p = n through- 
out the simulations, according to the procedure introduced in 
Section 2.3. 

In order to evaluate in quantitative terms the behaviour 
of the nested-grid technique, the model outlined above was 
executed with a two-level grid system (i.e., in nested-grid 
mode) as well as with a single-level grid (i.e., in single-grid 
mode) . The first grid level of the two-level grid system (hence- 
forth 2G), which covers the whole computational domain, had 
N r x = 147 x 455 grid zones (Ar ~ 1.46 x 10~ 2 and A(f> ~ 
1.40xl0~ 2 ), whereas the second level had N r x JV^ = 264x464 
zones (Ar ~ 7.3 x 10" 3 and A<j!> ~ 7.0 x 10~ 3 ). With this 
setup, the higher resolution region extends from r ~ 0.5 to 
r ~ 2.4 and from (j> — t/2 to (f> ~ 3ir/2. In order to achieve 
the same resolution with a single-level grid (henceforth 1G), 
the mesh must have N r x N$ — 290 x 904 grid zones. Although 
the grid system 2G covers nearly a half of the whole domain 
with the same numerical resolution as the grid 1G, the per- 
turbations induced by the planet propagate to the entire disc 
over a short time-scale and after 5 orbits the spiral wave pat- 
tern has already developed. Therefore, after a few orbits, one 
should expect that the results of the two simulations start to 
differ. The discrepancy depends upon the ability of the first 
level of the grid 2G to capture the same flow features as the 
grid 1G does, in the region (f> < tv/2 and <j> > 3 7r/2. 

Since this study is about migration, we focused on the 
evaluation and comparison of the semi-major axis evolutions, 
which also give a direct indication of the acting torques as 
a function of time. Note that this is a more strict test than 
simply comparing the torque distributions at certain times, 
which would only imply that a is the same at those times. 
The orbital decay for the two simulations is shown in Fig- 
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Fig. 15. — Semi-major axis evolution of the test model as 
calculated in a nested-grid mode with the grid system 2G 
(solid line) and in a single-grid mode with the grid 1G (dashed 
line). Radial and azimuthal resolutions (Ar ~ A(f> ~ 7.0 x 
10~ 3 ) coincide over roughly a half of the whole computational 
domain ([0.4,2.5] x 2vr). 
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Fig. 16. — Top. Normalised difference of the two functions 
02G and oig, shown in Figure 15, according to equation (Al). 
Each filled circle represents the average of the data over time- 
intervals of a tenth of an orbit that are centred at the middle 
point of each time-interval. Bottom. Running time average 
of the data displayed in the top panel, according to equa- 
tion (A2). 
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ure 15. The solid and dashed lines pertain to grids 2G and 
1G, respectively. To estimate the differences in more detail, 
we computed the normalised difference 



Aa 



= 2 



an - aj 
an + a\ 



(Al) 



where the labels I and II identify the used grids 1G and 2G, 
respectively. Since the time-step is different in the two sim- 
ulations, aiG and <i2G were averaged over time-intervals of 
0.1 orbits and the value measured from equation (Al) was 
assigned to the central time of each interval. As shown in 
the top panel of Figure 16, Aa/a is typically a few times 
10~ 5 and it doesn't increase beyond 1.5 x 10~ 4 . The bottom 
panel of Figure 16 illustrates the running time average of the 
normalised difference 



Aa 

a 



Aa 



(A2) 



which indicates that the discrepancies in the orbital evolu- 
tions are on average around a few times 10 -3 per cent. 

From the viewpoint of the computational load, the ad- 
vantage of the nested-grid technique is remarkable: in the 
computations reported above, the time-step required for nu- 
merical stability by grid 2G (first level) is twice as long as 
that required by grid 1G and so the time needed to complete 
one orbital period is twice as short. It is also worthwhile to 
point out that the refinement capabilities of the nested-grid 
strategy does not reduce the accuracy of the numerical algo- 
rithm, which strictly remains second-order accurate in space 
since the mesh step size is always constant on each grid level. 

For the sake of completeness, we show a test on the ac- 
celerated grid technique that we implemented and employed 
in this work. Other tests (not reported here) on the angu- 
lar momentum conservation of both the disc and the planet 
proved that conservation was achieved down to the machine 
precision. Thus, we refer to a more relevant situation and 
show how the migration calculated in a reference frame ro- 
tating with a variable Q compares to that calculated in a 
uniformly rotating grid with fi = fi = \J G M*/a^. The 
outcome of such comparison is illustrated in Figure 17. We 
computed the normalised difference (top panel) and the run- 
ning time average (bottom panel), where the quantity an in 
equation (Al) was obtained from the grid 1G with Q. =^ 
while the quantity a\ was obtained from the grid 1G with 
static rotation (i.e., Cl = 0). As Figure 17 proves, the dis- 
crepancy between the two semi-major axis evolutions is not 
significant. It has to be stressed that although the planet 
moves on a "continuous" path, the gravitational potential 
is centred in a grid zone. Therefore, results cannot be ex- 
actly the same since the planet's trajectory through the grid 
centres is different in the two simulations. In fact, in the 
model with Q, ^ there is only radial motion due to mi- 
gration, thus the time taken by the planet to cross a grid 
zone is Ar/a ~ (Ar/a) tm- Instead, in the other model, the 
azimuthal drift soon becomes the fastest component of the 

planet's motion and the time needed to cross a grid zone is 

-l 



then given by Acj>/ |fi K - fio| or (A0/27r) 1 - J (a/a f 
orbits. For instance, when Ar/a fa A<f> and a — 0.99 do, in 
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Fig. 17. — Comparison between the orbital evolutions com- 
puted with the grid 1G setting a constant rotation rate 
= Slo = \J G M*/ag and a variable rate (i.e., £7^0) such 
that the planet's azimuthal position remains constant. Top. 
Normalised difference (see eq. [Al]) sampled as in the top 
panel of Figure 16. Bottom. Running time average of the 
data displayed in the top panel. 



a steadily rotating grid a planet undergoes « tm/H ~ 300 
(see Fig. 15) as many encounters with the grid centres as it 
does in a grid rotating at a variable rate where there is no 
azimuthal drift. 



